{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "import datetime, os\n",
    "# hide tf logs\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'  # or any {'0', '1', '2'},\n",
    "# 0 (default) shows all, 1 to filter out INFO logs, 2 to additionally filter out WARNING logs, and 3 to additionally filter out ERROR logs\n",
    "import scipy.optimize\n",
    "import scipy.io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker\n",
    "import time\n",
    "from pyDOE import lhs         #Latin Hypercube Sampling\n",
    "import pandas as pd\n",
    "import seaborn as sns \n",
    "import codecs, json\n",
    "\n",
    "# generates same random numbers each time\n",
    "np.random.seed(1234)\n",
    "tf.random.set_seed(1234)\n",
    "\n",
    "print(\"TensorFlow version: {}\".format(tf.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Prep\n",
    "\n",
    "Training and Testing data is prepared from the solution file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_1 = np.linspace(-1,1,256)  # 256 points between -1 and 1 [256x1]\n",
    "x_2 = np.linspace(1,-1,256)  # 256 points between 1 and -1 [256x1]\n",
    "\n",
    "X, Y = np.meshgrid(x_1,x_2) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Data\n",
    "\n",
    "We prepare the test data to compare against the solution produced by the PINN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_u_test = np.hstack((X.flatten(order='F')[:,None], Y.flatten(order='F')[:,None]))\n",
    "\n",
    "# Domain bounds\n",
    "lb = np.array([-1, -1]) #lower bound\n",
    "ub = np.array([1, 1])  #upper bound\n",
    "\n",
    "a_1 = 1 \n",
    "a_2 = 4\n",
    "\n",
    "usol = np.sin(a_1 * np.pi * X) * np.sin(a_2 * np.pi * Y) #solution chosen for convinience  \n",
    "\n",
    "u = usol.flatten('F')[:,None] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot and save true solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.pcolor(x_1,x_2,usol, cmap = 'jet')\n",
    "# plt.axis('scaled')\n",
    "# plt.colorbar()\n",
    "# fig.savefig('Helmholtz_non_stiff_true.png', bbox_inches='tight', dpi = 500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trainingdata(N_u,N_f):\n",
    "    \n",
    "    leftedge_x = np.hstack((X[:,0][:,None], Y[:,0][:,None]))\n",
    "    leftedge_u = usol[:,0][:,None]\n",
    "    \n",
    "    rightedge_x = np.hstack((X[:,-1][:,None], Y[:,-1][:,None]))\n",
    "    rightedge_u = usol[:,-1][:,None]\n",
    "    \n",
    "    topedge_x = np.hstack((X[0,:][:,None], Y[0,:][:,None]))\n",
    "    topedge_u = usol[0,:][:,None]\n",
    "    \n",
    "    bottomedge_x = np.hstack((X[-1,:][:,None], Y[-1,:][:,None]))\n",
    "    bottomedge_u = usol[-1,:][:,None]\n",
    "    \n",
    "    all_X_u_train = np.vstack([leftedge_x, rightedge_x, bottomedge_x, topedge_x])\n",
    "    all_u_train = np.vstack([leftedge_u, rightedge_u, bottomedge_u, topedge_u])  \n",
    "     \n",
    "    #choose random N_u points for training\n",
    "    idx = np.random.choice(all_X_u_train.shape[0], N_u, replace=False) \n",
    "    \n",
    "    X_u_train = all_X_u_train[idx[0:N_u], :] #choose indices from  set 'idx' (x,t)\n",
    "    u_train = all_u_train[idx[0:N_u],:]      #choose corresponding u\n",
    "    \n",
    "    '''Collocation Points'''\n",
    "\n",
    "    # Latin Hypercube sampling for collocation points \n",
    "    # N_f sets of tuples(x,t)\n",
    "    X_f = lb + (ub-lb)*lhs(2,N_f) \n",
    "    X_f_train = np.vstack((X_f, X_u_train)) # append training points to collocation points \n",
    "    \n",
    "    return X_f_train, X_u_train, u_train "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PINN \n",
    "\n",
    "$W \\in \\mathcal{R}^{n_{l-1}\\times{n_l}}$ \n",
    "\n",
    "Creating sequential layers using the $\\textit{class}$ tf.Module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Sequentialmodel(tf.Module): \n",
    "    def __init__(self, layers, name=None):\n",
    "\n",
    "        self.W = []  # Weights and biases\n",
    "        self.parameters = 0 # total number of parameters\n",
    "\n",
    "        for i in range(len(layers)-1):\n",
    "\n",
    "            input_dim = layers[i]\n",
    "            output_dim = layers[i+1]\n",
    "\n",
    "            #Xavier standard deviation \n",
    "            std_dv = np.sqrt((2.0/(input_dim + output_dim)))\n",
    "\n",
    "            #weights = normal distribution * Xavier standard deviation + 0\n",
    "            w = tf.random.normal([input_dim, output_dim], dtype = 'float64') * std_dv\n",
    "\n",
    "            w = tf.Variable(w, trainable=True, name = 'w' + str(i+1))\n",
    "\n",
    "            b = tf.Variable(tf.cast(tf.zeros([output_dim]), dtype = 'float64'), trainable = True, name = 'b' + str(i+1))\n",
    "\n",
    "            self.W.append(w)\n",
    "            self.W.append(b)\n",
    "\n",
    "            self.parameters +=  input_dim * output_dim + output_dim\n",
    "            \n",
    "            \n",
    "    def evaluate(self,x):\n",
    "        \n",
    "        # pre-processing input \n",
    "        x = (x - lb)/(ub - lb) #feature scaling\n",
    "        \n",
    "        a = x\n",
    "        \n",
    "        for i in range(len(layers)-2):\n",
    "            \n",
    "            z = tf.add(tf.matmul(a, self.W[2*i]), self.W[2*i+1])\n",
    "            a = tf.nn.tanh(z)\n",
    "            \n",
    "        a = tf.add(tf.matmul(a, self.W[-2]), self.W[-1]) # For regression, no activation to last layer\n",
    "        \n",
    "        return a\n",
    "    \n",
    "    def get_weights(self):\n",
    "\n",
    "        parameters_1d = []  # [.... W_i,b_i.....  ] 1d array\n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "            \n",
    "            w_1d = tf.reshape(self.W[2*i],[-1])   #flatten weights \n",
    "            b_1d = tf.reshape(self.W[2*i+1],[-1]) #flatten biases\n",
    "            \n",
    "            parameters_1d = tf.concat([parameters_1d, w_1d], 0) #concat weights \n",
    "            parameters_1d = tf.concat([parameters_1d, b_1d], 0) #concat biases\n",
    "        \n",
    "        return parameters_1d\n",
    "        \n",
    "    def set_weights(self,parameters):\n",
    "                \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            shape_w = tf.shape(self.W[2*i]).numpy() # shape of the weight tensor\n",
    "            size_w = tf.size(self.W[2*i]).numpy() #size of the weight tensor \n",
    "            \n",
    "            shape_b = tf.shape(self.W[2*i+1]).numpy() # shape of the bias tensor\n",
    "            size_b = tf.size(self.W[2*i+1]).numpy() #size of the bias tensor \n",
    "                        \n",
    "            pick_w = parameters[0:size_w] #pick the weights \n",
    "            self.W[2*i].assign(tf.reshape(pick_w,shape_w)) # assign  \n",
    "            parameters = np.delete(parameters,np.arange(size_w),0) #delete \n",
    "            \n",
    "            pick_b = parameters[0:size_b] #pick the biases \n",
    "            self.W[2*i+1].assign(tf.reshape(pick_b,shape_b)) # assign \n",
    "            parameters = np.delete(parameters,np.arange(size_b),0) #delete \n",
    "\n",
    "            \n",
    "    def loss_BC(self,x,y):\n",
    "\n",
    "        loss_u = tf.reduce_mean(tf.square(y-self.evaluate(x)))\n",
    "        return loss_u\n",
    "\n",
    "    def loss_PDE(self, x_to_train_f):\n",
    "    \n",
    "        g = tf.Variable(x_to_train_f, dtype = 'float64', trainable = False)\n",
    "\n",
    "        k = 1    \n",
    "\n",
    "        x_1_f = g[:,0:1]\n",
    "        x_2_f = g[:,1:2]\n",
    "\n",
    "        with tf.GradientTape(persistent=True) as tape:\n",
    "\n",
    "            tape.watch(x_1_f)\n",
    "            tape.watch(x_2_f)\n",
    "\n",
    "            g = tf.stack([x_1_f[:,0], x_2_f[:,0]], axis=1)\n",
    "\n",
    "            u = self.evaluate(g)\n",
    "            u_x_1 = tape.gradient(u,x_1_f)\n",
    "            u_x_2 = tape.gradient(u,x_2_f)\n",
    "\n",
    "        u_xx_1 = tape.gradient(u_x_1,x_1_f)\n",
    "        u_xx_2 = tape.gradient(u_x_2,x_2_f)\n",
    "\n",
    "        del tape\n",
    "\n",
    "        q = -( (a_1*np.pi)**2 + (a_2*np.pi)**2 - k**2 ) * np.sin(a_1*np.pi*x_1_f) * np.sin(a_2*np.pi*x_2_f)\n",
    "\n",
    "        f = u_xx_1 + u_xx_2 + k**2 * u - q #residual\n",
    "\n",
    "        loss_f = tf.reduce_mean(tf.square(f))\n",
    "\n",
    "        return loss_f, f\n",
    "    \n",
    "    def loss(self,x,y,g):\n",
    "\n",
    "        loss_u = self.loss_BC(x,y)\n",
    "        loss_f, f = self.loss_PDE(g)\n",
    "\n",
    "        loss = loss_u + loss_f\n",
    "\n",
    "        return loss, loss_u, loss_f \n",
    "    \n",
    "    def optimizerfunc(self,parameters):\n",
    "        \n",
    "        self.set_weights(parameters)\n",
    "       \n",
    "        with tf.GradientTape() as tape:\n",
    "            tape.watch(self.trainable_variables)\n",
    "            \n",
    "            loss_val, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "            \n",
    "        grads = tape.gradient(loss_val,self.trainable_variables)\n",
    "                \n",
    "        del tape\n",
    "        \n",
    "        grads_1d = [ ] #store 1d grads \n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            grads_w_1d = tf.reshape(grads[2*i],[-1]) #flatten weights \n",
    "            grads_b_1d = tf.reshape(grads[2*i+1],[-1]) #flatten biases\n",
    "\n",
    "            grads_1d = tf.concat([grads_1d, grads_w_1d], 0) #concat grad_weights \n",
    "            grads_1d = tf.concat([grads_1d, grads_b_1d], 0) #concat grad_biases\n",
    "        \n",
    "        return loss_val.numpy(), grads_1d.numpy()\n",
    "    \n",
    "    def optimizer_callback(self,parameters):\n",
    "                \n",
    "        loss_value, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "        \n",
    "        u_pred = self.evaluate(X_u_test)\n",
    "        error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)\n",
    "        \n",
    "        tf.print(loss_value, loss_u, loss_f, error_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "N_u = 400 #Total number of data points for 'u'\n",
    "N_f = 10000 #Total number of collocation points \n",
    "\n",
    "# Training data\n",
    "X_f_train, X_u_train, u_train = trainingdata(N_u,N_f)\n",
    "\n",
    "layers = np.array([2, 3, 1]) #1 hidden layer\n",
    "\n",
    "PINN = Sequentialmodel(layers)\n",
    "\n",
    "init_params = PINN.get_weights().numpy()\n",
    "\n",
    "start_time = time.time() \n",
    "\n",
    "# train the model with Scipy L-BFGS optimizer\n",
    "results = scipy.optimize.minimize(fun = PINN.optimizerfunc, \n",
    "                                  x0 = init_params, \n",
    "                                  args=(), \n",
    "                                  method='L-BFGS-B', \n",
    "                                  jac= True,        # If jac is True, fun is assumed to return the gradient along with the objective function\n",
    "                                  callback = PINN.optimizer_callback, \n",
    "                                  options = {'disp': None,\n",
    "                                            'maxcor': 200, \n",
    "                                            'ftol': 1 * np.finfo(float).eps,  #The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol\n",
    "                                            'gtol': 5e-10, \n",
    "                                            'maxfun':  50000, \n",
    "                                            'maxiter': 10000,\n",
    "                                            'iprint': -1,   # no iteration updates\n",
    "                                            'maxls': 50})\n",
    "\n",
    "elapsed = time.time() - start_time                \n",
    "print('Training time: %.2f' % (elapsed))\n",
    "\n",
    "print(results)\n",
    "\n",
    "PINN.set_weights(results.x)\n",
    "\n",
    "''' Model Accuracy ''' \n",
    "u_pred = PINN.evaluate(X_u_test)\n",
    "\n",
    "error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)        # Relative L2 Norm of the error (Vector)\n",
    "print('Test Error: %.5f'  % (error_vec))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_u = 400 #Total number of data points for 'u'\n",
    "N_f = 10000 #Total number of collocation points \n",
    "\n",
    "# Training data\n",
    "X_f_train, X_u_train, u_train = trainingdata(N_u,N_f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building the Hessian matrix\n",
    "\n",
    "[Refer TensorFlow 2.3 Documentation Example](https://www.tensorflow.org/guide/advanced_autodiff#example_hessian)\n",
    "\n",
    "What is the 'block' matrix?\n",
    "\n",
    "A matrix of tuples representing the dimensions of sub-hessians and to track indices while assembling the Hessian using sub-hessians\n",
    "\n",
    "Sub-hessians: $ \\frac{1}{\\partial^{W^{[l]}}} \\big( \\frac{\\partial loss}{\\partial^{W^{[l]}}} \\big )$\n",
    "\n",
    "Each column contains the matrix $\\frac{\\partial loss}{\\partial^{W^{[l]}}}$\n",
    "\n",
    "Each row contains kernel $W^{[l]}$\n",
    "\n",
    "Each matrix element is a tensor of size (kernel[i].shape, kernel[j].shape) \n",
    "\n",
    "Example: kernel shape:(5,4), shape of matrix $\\frac{\\partial loss}{\\partial^{W^{[l]}}}$: (2,3)\n",
    "\n",
    "shape of sub-Hessian: (5,4,2,3)  \n",
    "\n",
    "We reduce the shape of theses higher order tensors into 2D tensors using tf.reshape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[(100, 100) (50, 100) (2500, 100) (50, 100) (50, 100) (1, 100)]\n",
      " [(50, 100) (50, 50) (2500, 50) (50, 50) (50, 50) (1, 50)]\n",
      " [(2500, 100) (2500, 50) (2500, 2500) (50, 2500) (50, 2500) (1, 2500)]\n",
      " [(50, 100) (50, 50) (50, 2500) (50, 50) (50, 50) (1, 50)]\n",
      " [(50, 100) (50, 50) (50, 2500) (50, 50) (50, 50) (1, 50)]\n",
      " [(1, 100) (1, 50) (1, 2500) (1, 50) (1, 50) (1, 1)]]\n"
     ]
    }
   ],
   "source": [
    "num_kernels = (len(layers)-1)*2 # total number of weight and bias tensors\n",
    "\n",
    "block = np.zeros((num_kernels,num_kernels),object) \n",
    "\n",
    "for j in range(num_kernels):\n",
    "    for i in range(j+1):\n",
    "        \n",
    "        if i == j:\n",
    "            s = tf.reduce_prod(PINN.W[i].shape)\n",
    "            block[i,j] = (s.numpy(),s.numpy())\n",
    "        else:\n",
    "            block[j,i] = (tf.reduce_prod(PINN.W[j].shape).numpy(), tf.reduce_prod(PINN.W[i].shape).numpy())\n",
    "            block[i,j] = block[j,i]     \n",
    "            \n",
    "print(block)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computation and assembly of sub-hessians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise Hessian \n",
    "# N x N square matrix , N = total number of parameters\n",
    "H_u = np.zeros((PINN.parameters,PINN.parameters)) \n",
    "\n",
    "# pointer to mark position of sub-hessian assembly \n",
    "pointer = np.array([0, 0]) # coordinates to upper left corner element of current block\n",
    "\n",
    "for j in range(num_kernels):\n",
    "    for i in range(j+1):\n",
    "\n",
    "        with tf.GradientTape() as tape2:\n",
    "            with tf.GradientTape() as tape1:\n",
    "                loss_value, loss_u, loss_f = PINN.loss(X_u_train, u_train, X_f_train)\n",
    "\n",
    "            g = tape1.gradient(loss_u, PINN.W[i]) #sub-gradient , n_in * n_out\n",
    "\n",
    "        h = tape2.jacobian(g,PINN.W[j]) # sub-hessian \n",
    "\n",
    "        #diagonal term\n",
    "        if i == j :\n",
    "            # reshape higher order tensor into 2D tensor\n",
    "            h_mat = tf.reshape(h, block[j,i]) # [?]\n",
    "            \n",
    "            # shape of block, block is square for diagonal terms\n",
    "            block_shape = h_mat.shape\n",
    "            \n",
    "            # Assemble block in H matrix\n",
    "            # position of assembly determined by 'pointer' and size of block\n",
    "            H_u[pointer[0]:pointer[0]+block_shape[0], pointer[1]:pointer[1]+block_shape[1]] = h_mat\n",
    "            \n",
    "            # move pointer to new poistion\n",
    "            # move to next row ---> determined by number of rows in current block\n",
    "            pointer[0] = pointer[0] + block_shape[0]\n",
    "            pointer[1] = 0\n",
    "            \n",
    "        #non-diagonal term    \n",
    "        else:\n",
    "            # reshape higher order tensor into 2D tensor\n",
    "            print(h.shape)\n",
    "            h_mat = tf.reshape(h, block[j,i]) \n",
    "            \n",
    "            # shape of block\n",
    "            block_shape = h_mat.shape\n",
    "            \n",
    "            # Assemble block in H matrix\n",
    "            # position of assembly determined by 'pointer' and size of block\n",
    "            H_u[pointer[0]:pointer[0]+block_shape[0], pointer[1]:pointer[1]+block_shape[1]] = h_mat\n",
    "            \n",
    "            # Assemble symmteric part by switching indices and transposing the block\n",
    "            H_u[pointer[1]:pointer[1]+block_shape[1], pointer[0]:pointer[0]+block_shape[0]] = tf.transpose(h_mat)\n",
    "            \n",
    "            # move pointer to new poistion\n",
    "            # move to next column ---> determined by number of columns in current block\n",
    "            pointer[1] = pointer[1] + block_shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise Hessian \n",
    "# N x N square matrix , N = total number of parameters\n",
    "H_f = np.zeros((PINN.parameters,PINN.parameters)) \n",
    "\n",
    "# pointer to mark position of sub-hessian assembly \n",
    "pointer = np.array([0, 0]) # coordinates to upper left corner element of current block\n",
    "\n",
    "for j in range(num_kernels):\n",
    "    for i in range(j+1):\n",
    "\n",
    "        with tf.GradientTape() as tape2:\n",
    "            with tf.GradientTape() as tape1:\n",
    "                loss_value, loss_u, loss_f = PINN.loss(X_u_train, u_train, X_f_train)\n",
    "\n",
    "            g = tape1.gradient(loss_f, PINN.W[i]) #sub-gradient , n_in * n_out\n",
    "\n",
    "        h = tape2.jacobian(g,PINN.W[j]) # sub-hessian \n",
    "\n",
    "        #diagonal term\n",
    "        if i == j :\n",
    "            # reshape higher order tensor into 2D tensor\n",
    "            h_mat = tf.reshape(h, block[j,i]) # [?]\n",
    "            \n",
    "            # shape of block, block is square for diagonal terms\n",
    "            block_shape = h_mat.shape\n",
    "            \n",
    "            # Assemble block in H matrix\n",
    "            # position of assembly determined by 'pointer' and size of block\n",
    "            H_f[pointer[0]:pointer[0]+block_shape[0], pointer[1]:pointer[1]+block_shape[1]] = h_mat\n",
    "            \n",
    "            # move pointer to new poistion\n",
    "            # move to next row ---> determined by number of rows in current block\n",
    "            pointer[0] = pointer[0] + block_shape[0]\n",
    "            pointer[1] = 0\n",
    "            \n",
    "        #non-diagonal term    \n",
    "        else:\n",
    "            # reshape higher order tensor into 2D tensor\n",
    "            h_mat = tf.reshape(h, block[j,i]) \n",
    "            \n",
    "            # shape of block\n",
    "            block_shape = h_mat.shape\n",
    "            \n",
    "            # Assemble block in H matrix\n",
    "            # position of assembly determined by 'pointer' and size of block\n",
    "            H_u[pointer[0]:pointer[0]+block_shape[0], pointer[1]:pointer[1]+block_shape[1]] = h_mat\n",
    "            \n",
    "            # Assemble symmteric part by switching indices and transposing the block\n",
    "            H_f[pointer[1]:pointer[1]+block_shape[1], pointer[0]:pointer[0]+block_shape[0]] = tf.transpose(h_mat)\n",
    "            \n",
    "            # move pointer to new poistion\n",
    "            # move to next column ---> determined by number of columns in current block\n",
    "            pointer[1] = pointer[1] + block_shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute Eigenvalues\n",
    "\n",
    "In this section we use the hermitian property of the Hessian and the approx Hessian inverse matrices and compute their eigenvalues using the 'tf.linalg.eigvalsh()' method and plot them for comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eigenvalues of Hessian (H) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_u = tf.convert_to_tensor(H_u, dtype = tf.float32)\n",
    "v_hess_u = tf.linalg.eigvalsh(H_u)\n",
    "\n",
    "H_f = tf.convert_to_tensor(H_f, dtype = tf.float32)\n",
    "v_hess_f = tf.linalg.eigvalsh(H_f)\n",
    "\n",
    "np.savetxt(\"hess_eigvals.txt\",  np.array([v_hess_u, v_hess_f]).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = np.loadtxt(\"prove_stiffness/split/hess_eigvals.txt\")\n",
    "s2 = np.loadtxt(\"prove_stiffness/split/hess_eigvals_stiff.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABFZUlEQVR4nO3de3zO9f/48cdrZ2aMTYiYcp4cZgjFHFJONSKWwscpUaqPSuf0DZ9UDgn1ITr9VHyQkKQyUU7ZHIYRNTnNYXPYZna4ruv1++N9WZttGNd1va9tz/vttpvrer1Pz+t9XZ7X+3q/X+/nS2mtEUIIUbp4mB2AEEII15PkL4QQpZAkfyGEKIUk+QshRCkkyV8IIUohSf5CCFEKSfIXQohSSJK/EEKUQl5mbFQpdQ8w0L79RlrrtmbEIYQQpZVy1B2+SqkFQE/gtNa6ca72+4H3AU/gY63127mmRQJVtNb/vdq6g4ODdUhIiEPiFEKI0iImJiZJa125oGmOPPL/FJgFfH65QSnlCcwG7gWOAb8rpVZorffZZ3kEGHatFYeEhLB9+3YHhiqEECWfUurvwqY57Jy/1noDcPaK5lbAIa31X1rrLOBr4EF7UDWBC1rr1ILWp5QaqZTarpTafubMGUeFKYQQAudf8K0OHM31/Ji9DYwj/k8KW1BrPVdrHa61Dq9cucBfLUIIIW6QKRd8AbTWb5i1bSGEKO2cnfyPA7flel7D3nbTsrOzOXbsGBkZGY5YnXADfn5+1KhRA29vb7NDEaLEc3by/x2oq5SqjZH0B2Bc5L1px44dIyAggJCQEJRSjlilMJHWmuTkZI4dO0bt2rXNDkeIEs9h5/yVUl8Bm4H6SqljSqlhWmsL8CTwAxAPLNZa73XE9jIyMggKCpLEX0IopQgKCpJfckK4iMOO/LXWUYW0rwZWO2o7uUniL1nk/RTCdUy74CuEEKJgFpuFb/d/y5/n/iSoTBDDwq55O1SRSfIXQgg3kmHJYMCSAXx74FsA6gfVd0ryl8JuLhQXF0fVqlWJi4szOxQhhBvKsmbRe1Fvvj3wLePbjQfg+bbPO2VbkvxdaPLkyWzatInJkyebHYoQws3YtI3hK4az5tAa5vWax9td3ubQU4foF9rPKdszq6qnB/AWUB7YrrX+zIw4XO2rr77K868QQgCsObSGf//wb+KT4pnYcSLDw4YDcEelO5y2TUd29VyglDqtlNpzRfv9SqkDSqlDSqkX7c0PYtzwlY1R8qFEWb58OSNGjKB///6sXbvW7HCEEG4sKT2Jvov7YtVWvn7oa16+52VWHljJQ4sf4u/zhdZlu2mOPO3zKXB/7oZcVT27AY2AKKVUI6A+sElr/W/gCQfG4BYiIyOZN28eH330EYsWLcppnz17Ns8888xNr//SpUt06NABq9UKGF0kx40blzP9vffeY8KECTe9HYBJkyYRGhpKkyZNaNasGVu3bgXg/PnzzJkzJ8+8bdsawzLMnDmThg0bMnDgwDzPo6KiaN++PRaLxSGxCVESTNs8jfTsdL7p/w39G/fnWMox/vXtv/jz7J9UKVfFads1q6rnMeCcfR5rQesrDlU9Q0NDqVatGnXq1Mn5CwgI4JVXXgFg4sSJjBkzJmf+3bt306RJk5ve7oIFC+jTpw+enp4A+Pr6smzZMpKSkm563blt3ryZVatWERsby+7du/npp5+47TajWkdByX/Tpk0AzJkzhx9//JGFCxfmef7VV1/RuXPnPF+IQpRmpy+e5oNtH9AvtB+NKjci25rNgKUDyLRmsrjfYvy8/Jy2bbOqei4D7lNKfQBsKGjB4lDVc/jw4Tz22GMcOnSIQ4cOcfDgQapWrcrQoUMZP3483bp1IywsLGf+uLg4hyT/hQsX8uCDD+Y89/LyYuTIkUyfPj3fvNOmTaNx48Y0btyYGTNmAHD48GEaNmzIiBEjCA0NpWvXrly6dCnfsomJiQQHB+Pr6wtAcHAwt956KwAvvvgif/75J82aNeP5543eCOXKlWPUqFH89ddfdOvWjenTp+d7HhkZmfOlIERptmTfEhrPaUymJZPX2r8GwCvrXmHT0U3M6zWPekH1nBuA1tphf0AIsCfX874Yo3ddfv4YMKuo623RooW+0r59+/K1uVpSUpKuWbOmzs7O1lprvW7dOt2lSxf9/vvv67CwMP3444/rDz/8UGuttc1m04GBgTo9Pf2mtpmZmamrVKmSp83f319fuHBB16pVS58/f16/++67+o033tDbt2/XjRs31mlpaTo1NVU3atRIx8bG6oSEBO3p6al37Nihtda6X79++osvvsi3rdTUVN20aVNdt25d/cQTT+j169fnTEtISNChoaH54tBa61q1aukzZ87ktOd+brFYdHBwcKGvzx3eVyGcLe5UnPb6Py8dPjdcx56I1VprbbFadKPZjfSolaMcth2MDjUF5tViW9Uzn58i8rfVfBjqjQZLOqzvnn/67UOMv4wk+LVv3mld1l9zk0FBQbRp04ZVq1YRGRnJ/PnzGT58OP3792fs2LF55k1ISKBKlSqUKVPmel9RgZKSkggMDMzXXr58eQYNGsTMmTNztvHrr7/Su3dv/P39AejTpw8bN27kgQceoHbt2jRr1gyAFi1acPjw4XzrLFeuHDExMWzcuJHo6Gj69+/P22+/zZAhQ244fk9PT3x8fEhNTSUgIOCG1yNEcWXTNkatGkUF3wqsfmQ1G49s5LYKtxFcNpilDy8lJDDEJXE4+7RPTlVPpZQPRlXPFU7epkuNGDGC+fPnc+HCBTZs2EDv3r0LnK+wUz7h4eE8+eSTdO7cmb179/LJJ5/w1FNP8eSTTzJ+/Ph885cpU6bQ4mfPPPMM8+fP5+LFi9eM+/KpHDAS8uWLsLNnz6ZZs2Y0a9aMEydO4OnpSUREBG+++SazZs1i6dKl11z3tWRmZuLn57xzmUK4o/TsdDp91gm/iX78dvQ3JkRMYPjK4Ty0+CE+jv0YgAbBDZx6nj83hx3526t6RgDBSqljwBta6/lKqctVPT2BBdpBVT3zudqRulfZq0/3C76uI/2CdOrUidGjRzN16lT69euHj48PAO+88w7JyclkZWUxffr0Ai/2Hj16lFatWjFr1iymT5/OkiVLuHDhAh988AEAWVlZ+bZXsWJFrFYrGRkZ+RJopUqVePjhh5k/fz5Dhw7lnnvuYciQIbz44otorfnmm2/44osvrvp6xowZk3OR+sCBAxw8eJC6desCsHPnTmrVqgVAQEAAqakFjsB5VcnJyQQHB0vNflHq/PuHfxN9OJqxrcYSXDaYd397l8S0RKZ1ncZTrZ9yeTzFuqqnO1BKMWTIEF599VX27DFucYiOjqZcuXK88MILDBgwAIDY2FhGjRqVZ9mYmBj++OMPRo0aRWJiIlWrVs3pKQTkfJFcqWvXrvz666906dIl37Rx48Yxa9YsAMLCwhgyZAitWrUCjAvUzZs3L/AUT0HS0tJ46qmnOH/+PF5eXtSpU4e5c+cCximvdu3a0bhxY7p168a77757XeuMjo6mR48e1zWvEMWd1pp5sfNY9ccqVv6xkufbPk/XO7oS+XUkFfwq8OvQX2lVvZV5wbn7n7te8L0sLS1Nx8XF5TwfM2aMHjx4sH788cd1q1at9IkTJ3StWrXyXex99dVX9c6dO7XWWvfp00f37dtXHz9+PGe6xWIpcHsxMTH60UcfdcIrcb7evXvrAwcOFDrdnd5XIW6GxWrRI1eM1ExA15lZRz/53ZM605KpD587rHss7KGPpxy/9kpuEiZe8C0V/P39ady4cc7zc+fOsXDhQo4cOUL9+vXp3r07c+bMyXexNyYmhqSkJDw8PHJugho3bhyVK1cmNTWV6dOnF3hxNywsjI4dO2K1WnP6+hcHWVlZREZGUq+ek7uwCWEyrTUjV45kwc4FvHz3y/Rp2IfZv8/G28ObWoG1WPXIKrNDlOTvDD179mTMmDH4+Phw6tQpypcvX+B8q1fnPxt2vXV/hg4delMxmsHHx4dBgwaZHYYQTmOxWdh6bCtfxn3Jgp0LeK39a9x3x310+rwTFXwrcCL1BNXLVzc7TMC8wm4RGIXd9gJfa63XmxGHs0RFRREVVeAlECFECXX20lkiv45k45GNADxQ7wGWxS/jrQ1vUadSHdYNWuc2iR/MK+ymgTTAjxJY2E0IUbqkZqZyzyf3sPX4Vj7s8SF/PvUn5zPPc8lyiXfvfZfNwzZzW4Xbrr0iF3Lkkf+nwCzg88sNuQq73YuR5H9XSq0ANmqtf1FKVQGmAQMdGIcQQrjU8z8+z/6k/fzw6A+0rt6aAN8AVj+ymtMXT1O7Ym2zwyuQKYXdtNY2+/RzgC9CCFFMrTiwgv/G/JenWz/N4r2LqT+rPhabBX8ff7dN/OD8c/4FFXZrrZTqA9wHBGL8WshHKTUSGAlQs2ZN50YphBBFdO7SOd7f+j5vbXiLO2+5k+jD0ew8uZNn73qWTEsmXj7u3Z/GlOi01sswKntebZ65wFyA8PBw7Yq4hBDierz1y1tM2jiJTGsmA+8ciKfy5IvdX7C8/3IebPDgtVfgBkpOYTchhHCypPQkZm6dyVsb3qJfo368ePeLpGamEvFZBC/d/VKxSfzg/OSfU9gNI+kPAB5x8jaFEMJhdp3cxWvRr7H9xHYS0xIBGNR0EJ88+AkeyoMMSwajw0czIWKCuYEWUckp7CaEEA50IvUE0zdPZ8bWGVT0q0j3ut1pUqUJzao2o0OtDkzbPI2hzYdSqUwlZveYbXa4RSaF3YQQws6mbTyz5hkW713M6YunUUrxaJNHmdZ1GkFlgwD4I/kPopZG8b99/8Pbw5un73ra5KhvjHtfji5h4uLiuPfee/nxxx+58847zQ5HCJHL6YuneXXdq8yLnUffRn1pWqUpUY2juKPSHQCkZaXxRvQbvL/1fXw8fZjQYQJjW4+9xlrdlyR/F5o8eTKbNm3ilVdeue4aPkII59p2fBujvxtNTGIMAC/d/RKTO0/ON99LP73ErN9nMSJsBG91fIsq5aq4OlSHUkbVT/cWHh6ut2/fnqctPj6ehg0bmhSRcBZ5X4UrJJxLYMIvE4hOiOZoylFuDbiVZ+96ltbVW3N3zbtRSmGxWdh5ciceyoOwamFcyLjAntN7aFezndnhXzelVIzWOrygaXLk7wTLly/nu+++IyUlhWHDhtG1a1ezQxJCYNyYNXnjZGZum4mn8iSyQSTNqzbn8fDHKe/7T/Xdb+K/YeiKoZzPOE/HkI6sG7yOCn4VilXivxbTkr9Syh/4BZigtTa/uLUDRUZGEhkZyblz53juuedykv/s2bM5ePAgM2bMMDdAIUqZS9mXmBc7jzd/eZNzl84xuNlgJnacWGCVzQ+2fsDTa56mZfWWPHvXs3Sq3cmEiJ3PkV09FwA9gdNa68a52u8H3sfo6vmx1vpt+6TxwGJHbd8MoaGhnD17Fn9//5y2U6dOMXbsWCZNmsTEiRNzxsMF2L17N61btzYjVCFKHa01m49t5rOdn7Fo7yIuZF6gc+3OvNf1PZpVbVbgMpuObmLsmrFENojkyz5fUsa7TIHzlQRmVfWsDuzDKOlcbA0fPpzExETeeecdwPiw1atXj6FDhzJ+/Hi6detGWFhYzvxxcXGMGDHCrHCFKDWW71/Om7+8yc6TOynrXZaHGj7E0OZD6VCrA0qpPPOmZ6fz4e8fMq7tONrUaMOXfb6kb6O+eHt6mxS9aziyn/8GpVTIFc05VT0BlFJfAw8C5QB/oBFwSSm1OlelT+zzXn9ht5hn4NzOm34NeVRsBi1mXHWWQYMGERYWxuTJk/Hy8mL9+vWEhITw3Xff8dNPP3HhwgUOHTrEqFGj0FoTHx9PaGioY+MUQuTxcezHjFg5ggbBDZjXax79Q/sT4BuQZx6rzcq+M/uIPhzNh9s/5EDSAYaHDaeCXwWi7iwdAzGZUtVTa/0kgFJqCJB0ZeKH4lHYLSgoiDZt2rBq1SoiIyOZP38+w4cPp3///owdm7f/b0JCAlWqVMk3jq8QwnH+t/d/jFw5kvvr3M+3A77Fx9Mn3zzRCdE88PUDpGWlAdAguAFrH1tLBb8Krg7XVKb29tFaf+qQFV3jCN2ZRowYwYwZM+jYsSMbNmxgwYIFBc4XFxdHkyZN8rW3bNmS8PBwDhw4wPLly+ncuTOtW7cmJSWFiIgIhg4dSosWLWjZsiUAw4YNy3kshPjHmkNrGLhsIG1va8vSh5fmS/yXsi9RxrsM7Wq2Y3T4aO6scidtb2vL7RVvNylic0lVz5vUqVMnRo8ezdSpU+nXrx8+PsYH7p133iE5OZmsrCymT5/O7t278yX/o0eP0r59e6ZOncqwYcM4d+4crVu3ZtYsY4iDjh070rlzZ1q2bMlHH33k8tcmhLvLtmZz+Pxh5u+Yz7TN0wi9JZRVj6yirHfZnHm01kzbPI23NrxFwtMJVCxTkSn3TjExavcgVT1vklKKIUOG8Oqrr7JnjzF8cXR0NOXKleOFF15gwIABAMTGxjJq1Kg8y8bExLB//36effZZ2rZty44dO2jRokXOdH9/f2JjY4mPj2fUqFFUrVqVCRMmuOy1CeHO1h9eT7//9SMpPQmAIc2G8N697xHoF5gzT2pmKk989wQL4xYS2SCSLGuWSdG6H6nq6QBjx46lV69eOXemLl26lLS0NHbv3k1CQgKJiYns2LGD9u3b51kuJiaGadOmUb9+fQBee+01+vXrB8CuXbuoWbMmsbGxzJgxg+bNm7v2RQnhhrTWfHvgW9b+uZaPYz+mTqU6vHvvu7S8tSWht+TtTPFV3Fe8Fv0aCecTmNhxIi/f83K+nj6lmVT1dAB/f38aN865tYFz586xcOFCjhw5Qv369enevTtz5szJd7F337591KlTJ+d5TEwMycnJ+Pr64unpyZQpU+jfvz+JiYl4eXlx33330bt3b5e9LiHcidaa5398nqmbp+Lv7U/vhr35b8//5jnSz21dwjq8Pb35Zcgv3F3zbtcGWwxIeQcn6NmzJ2PGjMHHx4dTp05Rvnz5AudbunRpnuerV+f/jiyoTYjSxqZtPP3908z6fRZPtnyS6fdPx8sjb/q6lH2JF396kfvr3E+3ut2Y3WN2gb19hEGSvxNERUURFVU6+goL4Wxaa55Y9QRzY+cyrs043r333Tynb06mneTTnZ8yf8d8Dp09RIBvAN3qdpPEfw2S/IUQbm3BjgXMjZ3Li+1eZHLnyXkS/9u/vs2bv7xJhiWDNjXaMKvbLO6rc5+J0RYfkvyFEG7rROoJxq0dR4daHZjUeVK+C7ahlUN5pPEjjL97PPWC6pkUZfFkSvJXSjUEngaCgZ+11h+aEYcQwj2lZKYw5dcpfLbrMzKtmczrNQ+AQ2cP8VXcV3h6ePLyPS/Tq34vetXvZXK0xZMpVT211vHAKKWUB0YhOEn+QggAthzbwoAlAziacpT769zP4KaDGbtmLL8d+Y3UrFQA+of2x6ZteCgPk6Mtvkyp6qm13qeUegB4AvjCgTEIIYqxS9mX6Pe/fnh5ePHb0N+4q8ZdvB79Ohv/3sjgpoNpXq05Laq1oFnVZtJn/yaZVdVzn9Z6BbBCKfUd8OWV6ytSVU8hRIkwc+tMjqUcY1ybcRxPMSrBTIiYwCN3PkKD4AYmR1eymFLVUykVAfQBfCnkBrDiUNVTCOE4W45tYdLGSfh4+DB181RCAkPo3bA3HspDEr8TmHLBV2u9HlhvxraFEO7Dpm0s3L2Qj2M/ZsORDXgoDwL9AvnxsR9pWqWpnNN3Imfv2RJf1bMo4uLiqFq1KnFxcWaHIoTpLmZd5Onvn2bQ8kEcTTlKgE8Azas2Z92gdYRVC8PTw9PsEEs0qerpQpMnT2bTpk288sorfPXVV2aHI4Qp0rPTmb1tNpN/ncz5jPM83Ohhvu77NRotR/ouJFU9XehywpfEL0qrL+O+5KnVT3E24ywAQWWCmNNjDkopFNJ7x5WkqqcTLF++nO+++46UlBSGDRtG165dzQ5JCNOdTDvJv5b/C6UUlcpU4qlWT/HInY8QVDbI7NDck80Kp38By0Wo4fgb2eQ3lhNERkYyb948PvroIxYtWpTTPnv2bJ555pmbXv+lS5fo0KEDVqsVgEmTJhEaGkqTJk1o1qwZW7du5fz588yZMyfPcm3bts15PHPmTBo2bMjAgQPzPM7KyqJ9+/ZYLJabjlOI3Mp6laVh5YZYbBZWDFjBhIgJUpKhMCfWwLe1YF1niJvglE1IbZ+bEBoaytmzZ/H3989pO3XqFGPHjmXSpElMnDiRMWPG5EzbvXs3rVu3vuntLliwgD59+uDp6cnmzZtZtWoVsbGx+Pr6kpSURFZWVk7yHz16dM5ymzZtynk8Z84cfvrpJ2rUqEGDBg1yHgN07tyZRYsWMXDgwJuOVQirzcrmY5t5ds2z7D61mw97fEi7mu3MDsv9WDMgaQscXwkHZkCFxuAfAnf/zznb01q7/V+LFi30lfbt25evzdWmTZumn3/++ZznNptN16lTRx86dEi/8MIL+scff8wzf5s2bfTvv/9+09tt06aNTkhI0FprvXTpUt2zZ8988/Tv31/7+fnppk2b6ueee05rrbW/v7/WWuvHH39ce3t768aNG2sfH5+cx9OmTdNaa71z507drVu3m47zRrjD+yocJ9uarVvNa6WZgPZ801Mv27fM7JDcS0aS1jvGa/1jhNZf+2m9EK0Xemi9sb/W2Wlan9msdcqhG149sF0XkldNSeZAJDAPWAR0vdb87pr8k5KSdM2aNXV2drbWWut169bpLl266Pfff1+HhYXpxx9/XH/44Ydaa+OLITAwUKenp9/UNjMzM3WVKlVynqempuqmTZvqunXr6ieeeEKvX79ea611QkKCDg0NzbPs5eSvtda1atXSZ86cyfdYa60tFosODg6+qThvlDu8r8Jx5m6fq5mAZgJ63A/jzA7HvVgtWv/UWesvPbVe00rrraO03thP60UVtI59ziGbuFryN6uw23JguVKqIvAesPZmtx/xaUS+todDH2Z0y9GkZ6fTfWH3fNOHNBvCkGZDSEpPou/ivnmmrR+y/prbDAoKok2bNqxatYrIyEjmz5/P8OHD6d+/P2PHjs0zb0JCAlWqVMk3lGNRJSUlERgYmPO8XLlyxMTEsHHjRqKjo+nfvz9vv/02ERERN7wNT09PfHx8SE1NJSAg4KbiFaXX8ZTjPPPDM3gqT+pUqsOkTpPMDsl9WC7Czpfh1M/Qej6UrQFbhsKlE3BbHwhx/ilX0wq72Wd51T692BoxYgQzZsygY8eObNiwgQULFhQ4X1xcHE2aNMnXHh4ezl133UV8fDwzZ85k27ZtxMbGorXG39+fKVOm5Jm/TJkyZGRk5Gnz9PQkIiKCiIgI7rzzTj777LObSv4AmZmZ+Pn53dQ6ROn20OKHSM9O57477uODbh/g6+VrdkjuIWkrrO8OWWehzkjj3P5P90BAXWj/DQS1dEkYphR2U0rFA28D32utYx2x/asdqZf1LnvV6cFlg6/rSL8gnTp1YvTo0UydOpV+/frh42MMHffOO++QnJxMVlYW06dPZ/fu3fmS/9GjR2nVqhWzZs1i+vTpLFmyhAsXLvDBBx8AkJWVlW97FStWxGq1kpGRgZ+fHwcOHMDDw4O6desCsHPnTmrVqkVAQACpqak39JqSk5MJDg7G29v7hpYX4mTaSXac3MFd1e9izaNrzA7HvcS9CR7ecO9vENwGUuLhlgho9yX4uq7bq7O7ehZU2K068BTQBeirlBpV0IJKqZFKqe1Kqe1nzpxxcpg3TinFkCFDmDRpEsOHDwcgOjqacuXKMWXKFBITEwGIjY2lZcu83+gxMTH88ccfjBo1ivXr13PixIk8XUEvf5FcqWvXrvz6668ApKWlMXjwYBo1akSTJk3Yt28fEyZMICgoiHbt2tG4cWOef/75Ir2m6OhoevToUaRlhLjsUvYl3trwFhabhS/6SMX2PFIPQeL3UPNhSNoEtmyo0Ag6/eDSxA/mFXabCcy8xjzFpqrn2LFj6dWrFw0bNgRg6dKlpKWlsXv3bhISEkhMTGTHjh20b98+z3IxMTFMnTqVpk2b8tBDD3H27Fm8vP55S6xWK56e+eubjBkzhunTp9OlSxdatGiRpwtnbl9+mbdSdlpaWs7jw4cPF/j48nJvv/32db12IQAyLZnMjZnLkvgl/HbkN6zaSt+GfalTqY7ZobkPy0X4fTSg4A/j1z3+taHmQ6aE4+zkXyoKu/n7+9O4cc41bs6dO8fChQs5cuQI9evXp3v37syZMyffxd6YmBiSkpLw8PCgYcOGREVFMW7cOCpXrkxqairTp0/Pc3H3srCwMDp27Fjol8PNyMrKIjIyknr15OYbcX2iE6IZvnI4f537i9DKoZT1LktZ77LM7l6sL+c51sUjsKoRWC+CV3m483Wo1d+40GsSKezmBD179mTMmDH4+Phw6tQpypcvX+B8q1fnr3pxvXV/hg4delMxFsbHx4dBgwY5Zd2i5Jm5dSb//uHf1KlUh7WPrmX7ie28vO5logdHc0u5W8wOzzzaBn9/DRcPQ+jL4BMM3gHg4Qnd94D/bddchbNJYTcniIqKIiqqwFJHQpQYi/Ys4uk1TxPZIJLPIz8nwDeAOdvnUC+oHi1ubWF2eOY4twsOfgin1kHqQQioD4lr/zm/f883bpH4QQq7CSGKKMOSwc9//czwlcNpe1tbFvddjLen0TNs+4nttK/V/hprKKGOrYSNDwEKfCuBdyCkHgBbJtR/FqrdC1W7mB1lDqntI4S4LlprPtn5CePWjuN8xnmqlavGor6LchK/TduY0GECIYEh5gbqajYrpP0Fx74BbQEUePrCLR2gei+o9bDRtdPNSPIXQlzT4fOHGblyJD/+9SPta7XnhbYvEBESgb/PP0UNPZQHw8KGmRili2kbbB0BR5dB9nlQnlB3FDSZaBz5uzlJ/kKIq9r490Z6fNkDjWZO9zk8Hv54gSNu7T61Gy8PLxpVbmRClC5izYCEzyH9BPz1KaT/DZ5loPHrUOdxKHur2RFeN0n+QohCrTm0hr6L+1KjfA3WPLrmqqd0Xot+jT+S/yB+TLzrAnS1pC2w7XHjsfI0Si73iAev4lcKxZTBXJRStyul5iullpixfSFE4eLPxPPm+jfpu7gv3RZ2IyQwhPVD1hea+NOy0vjPxv/wy+FfCL813LXBukJaAhyw35NaJQK67YZbe4LygojvimXiB8d29SxKVc+/gGGS/IVwL+cunaPLF104kXqCin4VefWeV3ml/Sv4XSXBvbbuNWZsnUHzqs0ZGTbShdE6mS0bjiyB7U8aR/l1RoKnH5zdBidWQdg0ozRDMWV2VU8hhBt59odnOZV2it9H/H7dR/E//vUjD9R/gG8HfOvk6FwkIwn2T4O/FkDGKQhsAvcsNRJ/ZjLEPANVOkL9p82O9KY47LSP1noDcPaK5pyqnlrrLOBr4MHrWV9xKexWFHFxcVStWpW4uDizQxEinwU7FvDZrs948e4Xi3T6ZvvI7czpPufaM7o7bfvn8YH3IagVdFgF98dCgL1GUeJasKRB0/9AARe9ixNTqnoqpYKUUh8BzZVSLxW0oNZ6rtY6XGsdXrlyZSeH6RqTJ09m06ZNTJ482exQhMiRYcng49iPGblyJF3v6MrrHV4v0vJ+Xn5UL1/dSdG5gM0K+96FdV2Nx37B0Ps4dFgB1XsYJRkuO7kWfCpBpeJ/bcOsqp7JQIGlnEuyy3V7rrd+jxDOtvf0Xjp+1pEz6WdoXb01S/otwcez4FLiV9Ja0/d/fXmg3gMMbjbYyZE6QXaKcSR/YAac+Q1q9AZrOngEgE9g/vm1Nuav2jnvF0IxJVU9nWD58uV89913pKSkMGzYMLp27Wp2SELkY9M2RqwcgU3b+OHRH+hcuzOeRUhqb//6Nsvil9G9Tv4hUt3eud3w491gSTXq6Lf5HEIeBaUKXyYl3hhmsWrJ+P8sVT2dIDIyksjISM6dO8dzzz2Xk/xnz57NwYMHmTFjhrkBilIvOT2Z6Vums/nYZj6L/IyudxQtoa39cy2vrHuFAY0HMLS5cyrMOtzJn+FCPNR/0uilc/sQqNkXgtuCx3WkwkT7UOPV7nVqmK4iVT1vQmhoKGfPnsXf/59b3E+dOsXYsWOZNGkSEydOZMyYMTnTdu/eTevWrc0IVYgc38R/w4ClA8iyZtGzXk8ea/JYkZa32qw8s+YZ6gXVY/4D81FXO1p2BxlJsO8/xkXcCqFQZ4RReyf8quNJ5XXpJOyfbizvX8t5sbqQVPW8CcOHDycxMZF33nkHMM6B1qtXj6FDhzJ+/Hi6detGWFhYzvxxcXGMGDHCrHCFYMuxLTyy7BGaV23OrO6zaFGtRZGT9/rD64lPimdx38WU9S7rpEgd4Mxm2PEcJG8DbFB7CLSYYST+orh4FDb2gaxkY4D1EkLKO9yEQYMGERYWxuTJk/Hy8mL9+vWEhITw3Xff8dNPP3HhwgUOHTrEqFGj0FoTHx9PaGio2WGLUigxNZG3NrzFvNh53Fb+NlZEreAW/xsbbKXz7Z3ZPmI7zas1d3CUDuZX2bg5q+HzEDIQAovwf89mhYsJcPAjODjH6AZ6zxKoFHbtZYuJEpH8n1nzDDtP7nToOptVbcaM+2dcdZ6goCDatGnDqlWriIyMZP78+QwfPpz+/fszduzYPPMmJCRQpUqVfEM5CuFMFzIu8O6md5m+ZTpZ1ixGhI3g9Q6v33Div8wtB2vRNjgVbYygFf6B0Tf/3g1FW0f6Mdj0GJxebzxXnlBrADSdVGJO91xWIpK/mUaMGMGMGTPo2LEjGzZsYMGCBQXOFxcXR5MmTfK1t2zZkvDwcA4cOMDy5cvp3LkzrVu3JiUlhYiICIYOHUqLFi1o2bIlAMOGDct5LMTV7E/az90L7ib5UjIDGg9gYseJ3FHpjptaZ3RCNB9s+4APun3gXn37L8TDxt6QcsAYRKXuE0U7Sr/4t3GU/+c8o3JnoxehTHW4LdLUcXadqUQk/2sdoTtTp06dGD16NFOnTqVfv374+Bh9pN955x2Sk5PJyspi+vTp7N69O1/yP3r0KO3bt2fq1KkMGzaMc+fO0bp1a2bNmgVAx44d6dy5My1btuSjjz5y+WsTxdvSfUtJvpTM1uFbaVW9lUPWueXYFr7Z/w2fRn7qkPU5xJnfYH0Po/xC2y/htt7G4+ulNUTfB6mHjJG2wmZAhQZOC9ddFO/7k92AUoohQ4YwadIkhg8fDkB0dDTlypVjypQpJCYmAhAbG5vviD0mJob9+/fz7LPP0rZtW3bs2EGLFv/8nPb39yc2Npb4+HhGjRrFhAkTXPa6RPG369Qubq94u8MSP8C+pH3UKF+D8r7lHbbOmxb3f+B3C9y3FUKiipb4wSjTnHIAWs2FjmtKReKHEnLkb7axY8fSq1cvGjZsCMDSpUtJS0tj9+7dJCQkkJiYyI4dO2jfPu/YpjExMUybNo369esD8Nprr9GvXz8Adu3aRc2aNYmNjWXGjBk0b+7mF9eE29l1ahdNqzR16Dr3nt5LaGU36bSQnQreAdB+uVGA7UbPyR/+whiQpWZfh4bn7iT5O4C/vz+NG+dUsebcuXMsXLiQI0eOUL9+fbp3786cOXPyXezdt28fderUyXkeExNDcnIyvr6+eHp6MmXKFPr3709iYiJeXl7cd9999O7d22WvSxRfF7MucjD5IAPvHOiwdVptVuKT4okIiXDYOm+IthmVNY+vgF6HwKsMlAu5sXVZM40LxDUiwduNfs24gCnJXynlD8wBsoD1WuuFZsThLD179mTMmDH4+Phw6tQpypcv+EO1dOnSPM9Xr85/O0RBbUJcS9zpODTaoUf+5zPO06RKE8KqmdTdMSMJEtfAseVwdCnUHWPU3L+eu3OvpLXx5bF3MmSdg9qDHB6uuzNlMBegD7BEa71SKbUIKFHJPyoqiqioAu95E8Ildp3cBUDTqo5L/kFlg9g6fKvD1lcklnRYexek/Wl0v2w6CRq9dPVaPIW5EA87x8PxlVCuDrT8CKrd5/iY3Zwpg7lgFHi7XNTe6sAYhBAY5/sr+FagVoUS0jfdq6xxs1bAHVD57qJf1AXjV8KWoXD4/4FnWWj+njEgy438cigBHFneYYNSKuSK5pzBXACUUpcHczmG8QWwE+lxJITD7Ti5g6ZVmzq07k73hd0p71uer/t+7bB1FspmhYTP4WyMcVReoxfUffzG1mXNhIyTsGM8HFkEjcZDg+eMuv2lmCmDuQDLgIeUUh8CKwtasCSO5CWEK6w8sJItx7bQpXYXh61z67GtfH/oe+6qcZfD1nlVfy2ArUPhr0/gyP9ufD3n42BFbfg2xEj8zd6BZm+X+sQP5g3mchH41zXmmQvMBQgPD9eFzOP+FQXFddO6wLdZXCeLzcKGvzfwr2//RbOqzXih3QsOW/d7m98j0C+QYc2HOWydhcpOg92vG6WW7/21aOf1s85B0ja4eNi4a/fPueDha/ThD6gLVSKcFXWxU2wHc/Hz8yM5OZmgoCD5AigBtNYkJyfj53cD53IFiamJdPi0AwfPHiTQL5CvHvoKX68iVq8sxOe7PmfpvqWMbzeeAN8Ah6zzqv6cb5ymuWfZ9Sf+9BOwbQQk/gDafhlReUHFptBukXGtQORRbAdzqVGjBseOHUNOCZUcfn5+1KhRMuuoOFNKZgo9vuzBidQTfNnnS3rW6+mwJJ1hyWDihol0vr0zr7R/xSHrLHhDScZRe/m6RnmGi4ehcpurL2OzwoU9RnmHvRONYRkbPm9cIwioA37VSsRwi85SbAdz8fb2pnbt2o5YlRDFjtaaM+ln+Pv83wz5dggHkg6wMmol3ep2c8j61xxaQ+vqralYpiLRg6Op7F/5usf2vS6nN8DZHZB5GtIS4OgyqNIROn4P/jWhxfR/5s06B8dXG3fxXjwM53dBZrJxWseSZsxTvgF0/AEC73RcjCWcKYO5FNL3XwhRCKvNyp7Te/jt6G/8dvQ3fv7rZ05dPAVA5bKVWfPoGrrcfnMXeG3axr4z+9h0dBNPrn6SwU0HM++BeddXvVNryEwyauiDcd79QpzRlnHGqI1vzYIIe/+OPRPh5I9Gn32fisaQiiGDIO2wMa5udooxXu7fi4z++LYsYznPslCxmZHsb4mA4DZQuS34h9xYn/9STLn6Ipu97/8f5Or7D0RprfcVtkx4eLjevn17kbd1Ke04WzeMIDkrk6TsDJKzM8myWbFqbf+zYUVj0xiPtcaKMc2Wa7/k3kM61zNdyDzGtH9a8y6fe57cy19tXdcRS75t2+fQ2j7VBlobyygfY5otE2wWwJZrPoX2tp8ysFxE27LzxqE8cm6D15Y0+/LGkkkWKxlake3hg8VmI9uaSRkPqOVjHGMcyMgmS6t/+lXbsinvqbjN2wsN7M/IworHPz/VbRYqeChu9Taex2dmo/EwEoZ9+YpeHtzi5YlGcyDDYsSXa3olTw+CvTywAocys41pysPYRzqbYC9PKnp6kK01h7MsYF+/RoO2UNnLkwoeiiyt+TvLal+/h/GKbVZu8fKgnKcHGTbN8Wz7dJR9f1mp4uWBv4fikk2TaLHa162M5bWVql6e+HkoLtpsnLHY0Hhi5fLn0YavUliBSzadc0OMn4cHt3h5E+xhYWBgGR4J9GNnhoVPz10C7/Jo5WV0b7RcZNatAVT28uCblEy+Op8BPhWMfWjNAGs6C6oHEODpwX/PXmLC6XROWmwANPLz5dfa/lT09LB/oGxGWYWgVsZrTPsL0o8a07QNtMV4TZcHN7+wFy5dvrznYfTL9y4HgU2NbWcmg/WicfNWdgrYMimQb2UIecQYjCWgnlHLR0nv8OullIrRWocXNM2M3j6F9f3Pk/yVUiOBkQA1a9a8oQ0lXUyk4+/f52tXGD85PBV4KoVHznNltAEeqDwHEoprP843Tf3Tcj3L51/2epc3Ere6nMA9vIx2bQVtyx+jh/3nu7agtC3/ej0zjMc2i7F8nv2gIMv4QkjKzuKiTRPibUyLz4JLGiAjZ32BHuCvjflPZUOWBpT9C0VDpg0CMJ6fsYBFW/NMt3lCoP35mWyAvNM9sVJRZaM1JFvyT/fDCvb8ddYCKFuuV6opr2x4K4VNa85bAWz8c9+hJsjDRlkPhdaQYrMnQWMvAJpbsVEB4zOUVsD0smiCPeCshnSbtsdnxKbQVFBWKnkozmjNUft0pcj5TN7vr6jupUjI1my6ZKOMpxdeeOChbKTZbPQsm01VrJzJtLDzkgUyLoBSxvuqLWRlpwGK05kW9mRYIOuCPT4r2CxYLemgFcpmIaKM4v7gmrQoF0B9j4t4Z52x7wpl/CkFWWeNL1Ctwaucvd0DI8H7QvZ5o82vilFp08ObPD3Ks1ONL4KAO4yDCK8AI6HnPC5vPPcKMH4RVGpuX4dwNDOO/PsC92uth9ufPwa01lo/WdgyN3rkn23NZuORjQSXDSaoTBBBZYPw9fQtOb2D/pgNsc8ady6C8Z+lfEOIWGU8PvE9nNtp/Cf18rf/Ww6qdzfmv3TSWPbytCKc0522eRrj1o4jIiSClVErKedTjgNJB/D08KSCbwUq+FVw7DliIUSRuduRv8t4e3rTqXYns8NwHGsmxE2A2x6CoHDjJ3SDf0P5RlC1M5S5Ne95z1u7GX+FKVP1hsI4lnKM16Nfp0fdHqyMWpnzZVo/uP4NrU8I4XpmJH+n9f0v8eLehH1vGz+Lg8LhlruNPxey2CyMWT0Gq7byQbcPSs6vKCFKGTOS/+9APaXUUuA80AXobkIcxUvaX7B/KoQ8BqEvmxbGmYtn2H5iO//p/B9qV5SutkIUVw5J/kUp56y1ttgTfxTGlcELjur7X2JZs+D3J407Fpv9x9RQqgVUI+6JOCqVqWRqHEKIm+OoPlOfAvfnbshVzrkb0AiIUko1sk9OBnpore8A9jgohpIr4XNI/B6aTYGy19Hn2knWHFpD3ClJ/EKUBA5J/lrrDcDZK5pzunRqrbOAy1064Z+Szg6LoUS7Yyh0+QXqF9ohyiVGrBzBpI2TTI1BCOEYzky8hZVzBinpfH1O/gQX9hn9qG9pf+35nehYyjGOpRyjTY1r1FsRQhQL13XOXyn1E1BQv8BXtNbfFnWjjirpXGIlbYGYZ+HsNuP29S4bTb91ffPRzQC0uU2SvxAlwXUlf631jRQNkS6dNyL1T/ilJ3j6Q+grUPcJ0xM/wOZjm/Hz8qNZ1WZmhyKEcABndvV0WjnnEsuaCb/0MmqldP7ZKEvrJrYc20KLai3krl0hSghHdfV0aTnnEsvT1zjSr9jcrRI/wJpH13D64mmzwxBCOIjLa/vciBut7VNsaJtR09wNRxv6YtcX1AqsRfta5l5wFkIU3dVq+0g3S3eQ8P9gVQNIdq8vuCxrFs/88AxzY+aaHYoQwsEk+ZvtbAzsHG+c6qkUZnY0efz010+cvXSWAY0HmB2KEMLBJPmbxXIRdr8Ba9sYg5u0nud2g1Qs2ruIQL9Aut7R1exQhBAO5l7ZpjQ5swn2/B/c1he67YKKTc2OKI9sazbL9y+nd4Pe0sNHiBKoRNfzdztaG6d5gsKh2r3QPQ4CG197ORMcPn8YPy8/etTtYXYoQggnkOTvKufjYOsIOLsdeuyD8vXcNvED1A2qy8lxJ7Fq67VnFkIUO5L8nS3zrFGHP/5d8A6EVvOgnPt16SyIUgovJR8RIUoiOefvTNYM+K4R7J0Mt/WDHnvhjn+Bh6fZkV3VpexLNJzdkMV7F5sdihDCSUw5rFNK3Q68AlTQWvc1IwansmYZg6F7+kHzd6FiMwi80+yorttvR39jf9J+AnwCzA5FCOEkDjnyV0otUEqdVkrtuaL9fqXUAaXUIaXUi5fb7TX+hzli224l5QDEvwdrwuDgf4222o8Vm8Svtebln1/mwa8fpJxPOe6u6drxgYUQrmPWSF4lizXLGFx99Z2w43nITnW72jzX40TqCf7z63+ICIlg09BNBPjKkb8QJZVDTvtorTcopUKuaM4ZyQtAKXV5JK9917NOpdRIYCRAzZo1HRGm82wbCQmfQa1HjNM8ZW81O6IbkmnNpH9of55v+zx3Vikev1aEEDfGlJG8lFJBSqmPgOZKqZcKWlhrPVdrHa61Dq9cubITw7wBWRdg12uQ8ofxvOHz0P5baLew2CZ+gNsr3s7Xfb+mxa0tzA5FCOFkZo3klQyMKupybuHE97B1GFxKBL/K9v76ocZfMZeWlUY5n3JmhyGEcAEZyet6nd0Bu1+FE6uhQmPjSD+opdlROVT43HBaVW/F570/NzsUIYSTOfO0T85IXkopH4yRvFY4cXvOdXC2MbZu00lw/+8lLvGnZ6dz8OxBagfWNjsUIYQLyEhe12LNMPrrt5oHzaeCTwWzI3KKPaf3YNM2mlZ1rwJzQgjncFRvn6hC2lcDqx2xDZezWYyqm4fmQvc94BdcYhM/wK6TuwBkgHYhSgkp3FKQtMPGRd1T6yBkIJTw4mZaa+bFzuMW/1sICQwxOxwhhAtI8s9N22DTo3BkEShvuOsTuH2I2VE5nVKKd+99l/TsdDzcbEAZIYRzSPIHyEwG3yBjJC3v8tDgOaj3JPjfdu1li7mzl85SqUwlOoR0MDsUIYQLmXaYp5SKVErNU0otUkq5fpzAi3/DvndhbVtYVuWfwdNbfQTNp5SKxH/4/GHqflCX/27/r9mhCCFczGHJ/waKuy3XWo/AuNmrv6PiuC6ZZ2F1U9j5gtGbp/EbUKb43pl7I7Kt2QxYMgCLzcK9d9xrdjhCCBdz5GmfT4FZQM4dQrmKu92LUd7hd6XUCq117vo+r9rncZ0jiyD7AnTZALfc49JNu4uXf36Zrce3srjvYm6veLvZ4QghXMxhR/5a6w3A2Suac4q7aa2zgMvF3VCGKcD3WutYR8VxXW7pAE0nQ+XSWbJ435l9vLf5PUa1GEW/0H5mhyOEMIGzL/gWVNyttf3xU0AXoIJSqo7W+qPcCzq1qmeFRsZfKRV/Jp4q/lX4v47/Z3YoQgiTXHfyd0Jxt5nAzKtMnwvMBQgPD9dFXX8+1iw48D6c3gANn4Mqpbd3y0ONHuLBBg/i5SGdvYQora77f3+xLu6W+if8NgDObofK7SAzyeUhuIvD5w9Tq0ItSfxClHLO7uppfnG3U9HwQytIPQT3LIV7f4WaD7k0BHeQbc3m1XWvcvv7tzN4+WCzwxFCmMxhh39uW9wteRuUqQodVkK50tmrJeFcAlFLo9h6fCtDmw3l/W7vmx2SEMJkSuubP53ubOHh4Xr79u03trDWYE0HL3/HBlVMxJ+Jp9PnnciwZDC351zp3SNEKaKUitFahxc0reSf+FWq1CZ+AIvNQnDZYBb1XUSjyqW3h5MQIq+Sn/xLmbSsNGZsmcGOkztY0m8Jd1a5k12jdknBNiFEHpL8S5D9Sfvpu7gv+87so3WN1qRkplDBr4IkfiFEPpL8S4BTaacYvXo0y/cvp6JfRdY+tpYut99Iz1whRGlhSvJXSjUEngaCgZ+11h+aEUdJUc6nHGlZaYxvN56xrcdStVxB9+IJIcQ/HDWG7wKgJ3Baa904V/v9wPsY3Tw/1lq/DaC1jgdGKaU8MArBSfIvgpTMFLYc28IPh35gXNtx3BpwK2sGrkEpZXZoQohiwlFH/p9SxIqeSqkHgCeALxwUQ6lw+uJpwueGczTlKB7Kg7BqYQxsMlASvxCiSBw1gPsGpVTIFc05FT0BlFKXK3rusy+zAlihlPoO+PLKdTq1sFsxZbFZiFoaxZn0M3zT/xvuqXkPQWWDzA5LCFEMOfOcf6EVPZVSEUAfwBdYXdDCDi/sVgKsPriadQnr+OTBT4hsEGl2OEKIYuy6kr8TKnquB9YXdbnSxKZtLN+/nJgTMaRmpTKz20x61evFz4N+plPtTmaHJ4Qo5q4r+Rfrip7FUPyZeIZ8O4Rtx7fhqTypFViLGffPwEN5SOIXQjiEM0/75FT0xEj6A4BHnLi9EmHfmX20W9AObw9vPn3wUwY2GSjll4UQDueQWz/tFT03A/WVUseUUsO01hbgckXPeGCxyyt6FiOXC+zVrVSXqMZRbBuxjcHNBkviF0I4Rcmv6unmtNbMi53HB9s+YN2gdVT2r2x2SEKIEuJqVT2l6IuJrDYro1aN4vFVj+Pv7c/pi6fNDkkIUUrIOQWTZFoyefSbR1mybwkv3/0yEztNlBu1hBAuI8nfJFN+m8KSfUuY1nUaz7Z51uxwhBCljCR/F7NpGx7Kg+fbPk9YtTB61utpdkhCiFLItHP+Sil/pdR2pVSpyH4ZlgwmrJ9A4zmNybBkUMa7jCR+IYRpHNXVc4FS6rRSas8V7fcrpQ4opQ4ppV68YrHxwGJHbN/d7UjcwV0f38Wbv7xJg+AGnM84b3ZIQohSzlFH/p8C9+duyFXVsxvQCIhSSjWyT7sXo8Bbie7ecj7jPN0WdiNsbhjHU4+zMmoly/ovk3r7QgjTmVXVMwLwx/hSuKSUWq21tuVeuCRU9fT28KZ6QHUmdpzI6JajqVimotkhCSEEYFJVT631KwBKqSFA0pWJ3z5Psa3qmW3N5vD5w9QNqsvHD3xsdjhCCJGPKVU9L9Naf3qjy7qjn//6mWXxy/jhzx+waRv7n9yPj6eP2WEJIUQ+UtXTAfae3sur0a+yfP9yAnwCaFm9Jc/e9SzeHt5mhyaEEAWSqp4OcOriKX766ycmdZrEuDbj8PXyNTskIYS4KkcN4P4VxkXcYKXUMeANrfV8pdTlqp6ewIKSVtVz2/FthN8aTvta7Ul+IVlO8Qghig1H9faJKqR9NYUM01jcbTu+jbbz2/JWx7d46Z6XzA5HCCGKRKp63oCk9CT6Lu5LjfI1GBU+yuxwhBCiyKS2TxGdTDvJQ4sf4tTFU/w29Dfpuy+EKJYk+RfB+Yzz3PXxXZy+eJoven9B+K0FjpEghBBuz5Tkr5SKAN4C9gJfa63XmxFHUQX6BfLzoJ9Jz07nzip3mh2OEELcMLMKu2kgDfDDuPPX7V3MugjAHZXukMQvhCj2TCnsBmzUWnfDqOz5poNicIpL2ZdYsGMBDWY3YNi3w8wORwghHMIhyV9rvQE4e0VzTmE3rXUWcLmwG7lq+ZwD3PaOqBlbZlBjeg2GrRhGoF8gA5sMNDskIYRwCFMKuyml+gD3AYHArIIWNruq55ELRxi3dhwdanXgjQ5v0L5WexljVwhRYphS2E1rvQxYdo15TK3q6aE8eCL8CV5o9wI1KxTPktJCCFEYKexWiBrlazCre4E/SoQQothz5h2+OYXdlFI+GIXdVjhxezfNpm28uu5V2i1ox7qEdWhdrIYREEKI6+aorp5fAZuB+kqpY0qpYVprC3C5sFs8sNjdC7u9t+k9Jm2chJeHF5uObsKWf4wZIYQoEaSwm93aP9fy8s8v069RPxb1XSQXd4UQJZoUdgNSM1N5ZOkjNKrciI8f+FgSvxCixCvVtX1SM1Mp612WAN8AVkStoH5Qfcr7ljc7LCGEcLpSe+S/6o9VNJrTiDm/zwGg7W1tCSobZHJUQgjhGqUu+SelJxG1NIpeX/Wigm8FWlZvaXZIQgjhcmZV9fTAqOpZHtiutf7MFdvNtmbz4NcP8vvx33kz4k1evPtFGXpRCFEqmVXV80GMm76ycWFVz9jEWGJOxPB57895vcPrkviFEKWWWVU96wObtNb/Bp5wUAzX1LpGaw6NPcSAxgNctUkhhHBLplT1xDjaP2d/bC1onUqpkUqp7Uqp7WfOnLmp+E6knuDrPV8DRtkGIYQo7Zx5wbegqp7V7Y+XAfcppT4ANhS0sNZ6rtY6XGsdXrly5ZsK5NV1rzL026Ekpibe1HqEEKKkMKuqZzrgkpFRMiwZLI1fSv/G/akWUM0VmxRCCLdX4qt6fn/we1IyU4hqXGAFCiGEKJVKfFXPr/Z8xS3+t9CpdidXb1oIIdxWia7qadM29p3ZR79G/fDyKNWVLIQQIo8SXdXTQ3kQ90Qc6dnpZoUghBBuqcSXd1BK4e/jb3YYQgjhVkp88hdCCJGfJH8hhCiFJPkLIUQpZFZVz3uAgfbtN9JatzUjDiGEKK1Mqeqptd6otR4FrAJcUs5ZCCHEP8yq6nnZI8CXDopBCCHEdTKrqidKqZrABa11akHrdGRVTyGEEHmZVdUTjMJunxS2sCOregohhMjLlKqeAFrrN6533piYmCSl1N83sh27YCDpJpZ3FomraCSuopG4iqYkxlWrsAnFoqqn1vqmDv2VUtu11uE3sw5nkLiKRuIqGomraEpbXCW+qqcQQoj8SnRVTyGEEAUr0VU9c5lrdgCFkLiKRuIqGomraEpVXEpr7Yz1CiGEcGNS20cIIUohSf5CCFEKlejkX1htIRPiuE0pFa2U2qeU2quUetrePkEpdVwptdP+192E2A4rpeLs299ub6uklPpRKXXQ/m9FF8dUP9c+2amUSlFKPWPW/iqodlVh+0gZZto/c7uVUmEujOldpdR++3a/UUoF2ttDlFKXcu23j5wR0zViK/S9U0q9ZN9fB5RS97k4rkW5YjqslNppb3fJPrtKbnD+50trXSL/AE/gT+B2wAfYhVFB1IxYqgFh9scBwB8Y9Y4mAM+ZvJ8OA8FXtL0DvGh//CIwxeT38STGzSqm7C+gPRAG7LnWPgK6A98DCrgL2OrCmLoCXvbHU3LFFJJ7PpP2V4Hvnf3/wS7AF6ht/z/r6aq4rpg+FXjdlfvsKrnB6Z+vknzkf9XaQq6ktU7UWsfaH6didH2tfvWlTPUg/1Rb/QyINC8UOgN/aq1v5g7vm6ILrl1V2D56EPhcG7YAgUqpaq6ISWu9VhtdrAG2YNxY6XKF7K/CPAh8rbXO1FonAIcw/u+6NC6llAIeBr5yxravElNhucHpn6+SnPyvVVvIFEqpEKA5sNXe9KT959sCV59esdPAWqVUjFJqpL2titY60f74JFDFhLguG0De/5Bm76/LCttH7vK5G4pxhHhZbaXUDqXUL8oYT8MMBb137rK/7gFOaa0P5mpz6T67Ijc4/fNVkpO/21FKlQOWAs9orVOAD4E7gGZAIsbPTle7W2sdhlF6e4xSqn3uidr4rWlKf2Bl3Bn+APA/e5M77K98zNxHBVFKvQJYgIX2pkSgpta6OfBv4EulVHkXh+WW710uUeQ9yHDpPisgN+Rw1uerJCd/h9UWcgSllDfGm7tQa70MQGt9Smtt1VrbgHk46efu1Witj9v/PQ18Y4/h1OWfkvZ/T7s6LrtuQKzW+pQ9RtP3Vy6F7SNTP3dKqSFAT2CgPWlgP6WSbH8cg3FevZ6rYrJvt7D3zvT/p0opL6APsOhymyv3WUG5ARd8vkpy8neb2kL284nzgXit9bRc7bnP1fUG9ly5rJPj8ldKBVx+jHHBcA/Gfhpsn20wcEOVWx0gz9GY2fvrCoXtoxXAIHuvjLswxqxILGgFjqaUuh94AXhAa52eq72yMgZXQil1O1AX+MsVMeWKobD3bgUwQCnlq5SqbY9tmytjA7oA+7XWxy43uGqfFZYbcMXny9lXs838w7gy/gfGt/YrJsZxN8bPtt3ATvtfd+ALIM7evgKo5uK4bsfoabEL2Ht5HwFBwM/AQeAnoJIJ+8wfSAYq5GozZX9hfAElAtkY51iHFbaPMHphzLZ/5uKAcBfGdAjjfPDlz9hH9nkfsr+/O4FYoJcJ+6vQ9w54xb6/DgDdXBmXvf1TYNQV87pkn10lNzj98yXlHYQQohQqyad9hBBCFEKSvxBClEKS/IUQohSS5C+EEKWQJH8hhCiFJPkLIUQpJMlfCCFKof8P9zh9s7U8HsYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "### Plotting ###\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "\n",
    "plt.plot(s1[:,0], '--', color = 'orange', label = r'$\\nabla_{\\theta}^2 \\hat J_{BC}$ (Non-Stiff)')\n",
    "plt.plot(s1[:,1], color = 'orange', label = r'$\\nabla_{\\theta}^2 \\hat J_{PDE}$ ')\n",
    "plt.plot(s2[:,0], '--', color = 'green', label = r'$\\nabla_{\\theta}^2 \\hat J_{BC}$ (Stiff)')\n",
    "plt.plot(s2[:,1], color = 'green', label = r'$\\nabla_{\\theta}^2 \\hat J_{PDE}$')\n",
    "plt.yscale('symlog')\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig('prove_stiffness/split/split_200.png', dpi = 500)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
